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We apply a recently proposed novel thermostating mechanism to an interacting many-particle 
system where the bulk particles are moving according to Hamiltonian dynamics. At the bound- 
aries the system is thermalized by deterministic and time-reversible scattering. We show how this 
scattering mechanism can be related to stochastic boundary conditions. We subsequently simulate 
nonequilibrium steady states associated to thermal conduction and shear flow for a hard disk fluid. 
The bulk behavior of the model is studied by comparing the transport coefficients obtained from 
■ computer simulations to theoretical results. Furthermore, thermodynamic entropy production and 

exponential phase-space contraction rates in the stationary nonequilibrium states are calculated 
showing that in general these quantities do not agree. 
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I. INTRODUCTION 



Driving macroscopic systems out of equilibrium requires external forces. Now, the very existence of a nonequilib- 

■ rium steady state implies that the temperature of the system must remain time independent. One way to prevent the 
system from heating up indefinitely in nonequilibrium is the introduction of a thermostating algorithm |l|] . Starting 

£f) • from molecular dynamics simulations Evans, Hoover, Nose and others proposed deterministic thermostats to model 
equilibrium and nonequilibrium fluids (^-0] . In this formalism the (average) internal energy of the dynamical system 

0^ , is kept constant by subjecting the particles to fictitious frictional forces, thus leading to microcanonical or canonical 
distributions in phase space |q-|l~l|]. The major feature of this mechanism is its deterministic and time-reversible char- 
acter, which is in contrast to stochastic thermostats lH,|l^-|To|l ■ This allows to elaborate on the connection between 
microscopic reversibility and macroscopic irreversibility and has led to interesting new links between statistical physics 
and dynamical systems theory |^,P, |ll|^7| , especially to relations between transport coefficients and Lyapunov expo- 
q ■ nents ]l8|-p2[, and between entropy production and phase-space contraction p3| , |2^ . pi)| , [r5| , p^U2^ p8[ . Although used 
almost exclusively in the context of nonequilibrium systems, the abovementioncd thermostating mechanism presents 
• the drawback that the dynamical equations themselves are altered, even in equilibrium. This raises the question of 
whether some of the results are due to the special nature of this thermostating formalism or are of general validity 

Recently, an alternative mechanism in which thcrmalization is achieved in a deterministic and time-reversible way 
has been put forward by Klages et al. |54],[35| and has been applied to a periodic Lorentz gas under an external 

■ field. In the present paper we apply this thermostating method to an interacting many-particle system subjected 
to nonequilibrium boundary conditions giving rise to thermal conduction and to shear flow. The model is closely 
related to that of Chernov and Lebowitz p5|Jl6|j3qj37[] , who study a hard disk fluid driven out of equilibrium into a 
steady state shear flow by applying special scattering rules at the boundaries in which the particle velocity is kept 
constant. Our model is introduced in Section || where the thermalization mechanism is tested under equilibrium 



conditions. In Section HI we move on to the case of an imposed temperature gradient and a velocity field by adapting 



the scattering rules, and we compute the respective transport coefficients. Having a deterministic and time-reversible 



system at hand we proceed in Section IV to investigate the relation between thermodynamic entropy production and 
phase-space contraction rate in nonequilibrium stationary states. The main conclusions are drawn in Section [y| 



II. EQUILIBRIUM STATE 

Consider a two-dimensional system of hard discs confined in a square box of length L with periodic boundary 
conditions along the x-axis, i.e., the left and right sides at x = ±L/2 are identified. At the top and bottom sides of 
the box, y = ±L/2, we introduce rigid walls where the discs are reflected according to certain rules to be defined later. 
The discs interact among themselves via impulsive hard collisions so that the bulk dynamics is purely conservative. 
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In the following and in all the numerical computations we use reduced units by setting the particle mass to, the disk 
diameter a and the Boltzmann constant ks equal to one. 

Before proceeding to the nonequilibrium case we define the disc- wall collision rules in equilibrium and check whether 
the system is well-behaved. Now, in equilibrium the bulk distribution is Gaussian with a temperature T, and the in- 
and outgoing fluxes at the top and bottom wall have the form (see |12-HlJ 




$(v x ,vy) = (2 7 rT 3 )- 1 /2|^| exp y-^fLj , (I) 

with v y < for the bottom wall and v y > for the top wall. Imposing stochastic boundary conditions on the systems 
in this setting would mean that for every incoming particle the outgoing velocities are chosen randomly according 
to Eq. (0). In practice, this is usually done by drawing numbers from two independent uniformly distributed 
random generators £, £ £ [0, 1] and then transforming these numbers with the invertible map 1 1 : [0, 1] x [0, 1] — > 

[0, oo ) x [0, oo) as 



(v x ,v y ) = T = VZT err 1 (CIV^W)) , (2) 



which amounts to transforming the uniform densities p(Q — 1 and p(£) = 1 onto <&(v x ,v y ) according to 



d(d£ 




dT (v x ,v y ) 


dv x dv y 




dv x dv y 



vl + vl 



^(2/ 7 rT 3 ) 1 /2 K |exp^^^j. (3) 

Note that so far we have restricted Eqs. (^|),(|^) to positive velocities v x ,v y € [0, oo), which implies a normalization 
factor in Eq. (||) being different to the one of Eq. (|j). In analogy to stochastic boundaries we now define the 
deterministic scattering at the walls as follows. First, take the incoming velocities v x , v y and transform them via 
T (v Xl v y ) = (C)£) onto the unit square. Second, use a two-dimensional, invertible, phase-space conserving chaotic 
map M : [0, 1] x [0, 1] — > [0, 1] x [0, 1] to obtain (£',£') = M(CiQ- Finally, transform back to the outgoing velocities 
via (v' x , v'y) = 1 1 (C', £,')■ I n order to render the collision process time-reversible, we also have to distinguish between 
particles with positive and negative tangential velocities by using Ai and M , respectively. Thus, particles going in 
with positive (negative) velocities have to go out with positive (negative) velocities and the full collision rules read 

(v' x ,v' y ) = T - 1 oMoT (v x ,v y ), v x >0 (4) 
(v' x ,v' y ) = T -'oM-'oT (v x ,v y ), v x <0, 

where T is meant to be applied to the modulus of the velocities |g . Since both the positive and the negative side of 
the tangential velocity distribution of Eq. (Q) are normalized to 1/2, this normalization factor has to be incorporated 
in Eq.([|) to render the full desired flux $ equivalent to the one of Eq. ([[]). Rewriting Eq. (||) in polar coordinates 
yields precisely the transformation used in [ p4|3^ |, in the limiting case where it mimics a reservoir with infinitely 
many degrees of freedom. It should also be realized that for obtaining the transformation T the total number of 
degrees of freedom of the reservoir has been projected out onto a single velocity variable, which couples the bulk to 



the reservoir. Eckmann et al. |39| used a similar idea to go from a Hamiltonian reservoir with infinitely many degrees 
of freedom to a reduced description when modeling heat transfer via a finite chain of nonlinear oscillators. It remains 
to assign the form of the chaotic map M. and we shall first adopt the choice of a baker map, as in Refs. p^j35|] , 

(C £>) - M(C £) - I (2C ' C/2) ; °^<l/2 f5l 

Later on we will investigate the consequences of choosing other mappings like the standard map (see, e.g., pfl| , ^l[ |). 
Since in equilibrium the in- and outgoing fluxes have the same form as in Eq. ([!]), the baker map yields a uniform 
density and our scattering prescription in Eq.(|]) can be viewed as a deterministic and time-reversible counterpart of 
stochastic boundary conditions. M. being chaotic, the initial and final momentum and energy of any single particle 
are certainly different, but both quantities should be conserved on the average. The latter is confirmed by numerical 
experiments in equilibrium, where, as usual in hard disk simulations, we follow a collision-to-collision approach [0. 
Keeping the volume fraction occupied by N — 100 hard disks equal to p = 0.1 sets the length of the box equal to 
L = 28.0. After some transient behavior which depends on the temperature of the initial configuration the bulk 
distribution is Gaussian with zero mean and mean kinetic energy T/2 in each directions. The in- and outgoing fluxes 
at the walls are correctly equipartioned with T as well and have the desired form of Eq. ([!]), so the system reproduces 



2 



the correct statistic properties. We close this section by a remark on how we measure the temperature of a flux 
to or from the boundaries. As the temperature of the tangential component we use the variance of the velocity 
distribution, T x := {{v x — {v x ) x ) 2 ) , where denotes an average over the density p(v x ). On the other hand, since 
in the normal direction we actually measure a flux, the appropriate prescription to measure the temperature of this 
component is T y := [v y ] y / [vy 1 ] , where [] y represents an average over the flux <I> and the denominator serves as a 

normalization. The temperatures of the in- and outgoing fluxes at the wall are then defined as Tii '■= (T x + T y )/2, 
tmdT w := (T l + T )/2. 

III. NONEQUILIBRIUM STEADY STATE 
A. Heat flow 

1. The Model 

In the following, we explicitly indicate the dependence of the transformation T on the parameter T by writing 
Tt- This immediately indicates how we may drive our system to thermal nonequilibrium: We just have to use 
different values of this parameter for the upper (T u ) and the lower wall {T d ). We deliberately avoid to use the word 
'temperature' for this parameter, since in contrast to stochastic boundary conditions we have generally no idea how a 
different T affects the actual temperature of the wall in the sense of the definition given above. In a nonequilibrium 
situation the temperature of the ingoing flux generally does not match exactly the parameter T. Therefore, we 
do not transform onto the uniform invariant density of the baker map anymore, and consequently, the outgoing flux 
might have all kinds of shapes or temperatures. Nevertheless, the hope is that the mapping M. will be chaotic enough 
to smooth out most of the differences and to produce a reasonable outgoing flux <E> such that the system is correctly 
thermostated. And this is indeed what we find in the numerical experiments. 

2. Numerical Results 

We set T u = 2, T d = 1 and p = 0.1 and average over about 40000 particle-particle collisions per particle and 
about 6000 particle-wall collisions per particle. We divide the available vertical height L — 1 into 20 equally spaced 
horizontal layers and calculate the time averages of the number density n(y) of the particles, the mean velocities 
u x (y) = (v x ), u y (y) — (v y ), and the variances ((v x — u^) 2 ), ((v y — u y ) 2 y Furthermore, we record the time average 
of the kinetic energy transfer and measure the temperatures of the in- and outgoing fluxes of both walls as described 
in the preceding section. The temperatures at the walls are then defined as the mean value of the in- and outgoing 
temperatures, T % J d := (T" /d + T u/<i )/2. Time series plots of these quantities confirm the existence of a nonequilibrium 
stationary state (NSS) induced by the temperature gradient. 

Fig. [I] shows the temperature profile between the upper and the lower wall. Apart from boundary effects it is 
approximately linear, and the respective kinetic energy is equipartitioned between the two degrees of freedom. The 
parameters T u and T d are represented as (*), and we find a 'temperature' jump, whereas the measured temperatures 
(+) at the walls seem to continue the bulk profile reasonably well. The profile of the number density n = 4/irp 
is depicted in Fig. ||. Note again the boundary effects. The densities of the in-coming particles at the upper wall are 
Gaussian shaped (Figs. ||a,c), whereas the outgoing densities (Figs. ||b,d) show cusps due to the folding property of 
the baker map. Nevertheless, the baker map produces a reasonable outgoing flux which generates a NSS. 

In order to examine the bulk behavior we now compute the thermal conductivity in our computer experiment and 
compare it to the theoretical value. For this purpose, we measure the heat flux Q across the boundaries and estimate 
the temperature gradient dT(y)/dy by a linear least square fit to the experimental profile. To discard boundary effects 
we use only data in the bulk of the system, namely from layer 3 to layer 18, i.e., excluding the top two and the bottom 
two layers. The experimental heat conductivity is then defined as 

•W = Q (Jj^J > ( 6 ) 

whereas the theoretical expression for the conductivity of a gas of hard disks with unit mass and unit diameter as 
predicted by Enskog's theory reads [|^|,fl3| 
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A, = 1.0292 



X 2 



6n + 0.8718(6n) 2 x 



(7) 



Here, b is the second virial coefficient, b = tt/2, and \ 1S the Enskog scaling factor, which is just the pair correlation 
function in contact [E3, 



X = 



(l-fn)2- 



(8) 



Since (0) depends on local values of T and n we define the theoretical effective conductivity Xth as the harmonic mean 
over the layers E2], 



iV,, 



At?i — I i/Ni ayers 2^ 



(9) 



Table || compares A eKp to Ath by showing the ratio of the experimental to the theoretical conductivity for different 
particle numbers and temperature differences. The agreement is quite good, so our thermostating mechanism produces 
a NSS which is in agreement with hydrodynamics. 

Furthermore, going into the hydrodynamic limit by increasing the number of particles we observe that the discon- 
tinuity in the outgoing flux of Figs. ||(b,d) diminishes, as expected, since both the in- and outgoing flux come closer 
to local equilibrium. 



B. Shear Flow 

Inspired by the recently proposed model of Chernov and Lebowitz [T^ , p^ ] for a boundary driven planar Couette- 
flow in a nonequilibrium steady state, we now proceed to check whether it is possible to combine our thermostating 
mechanism with a positive (negative) drift imposed onto the upper (lower) wall, respectively. Chernov and Lebowitz 
chose a purely Hamiltonian bulk and simulated the drift at the boundaries by rotating the angle of the particle velocity 
at the moment of the scattering event with the wall while keeping the absolute value of the velocity constant. This 
setting could be formulated in a time-reversible way and keeps the total energy of the system generically constant. 
Here we separate the thermostating mechanism and the drift of the walls by introducing the map 

Sd(v x ,v y ) = (v x +d,v y ), (10) 

and by applying this shift to the 'thermostated' velocities. Time-reversibility forces us to do the same before thermo- 
stating. Thus, the full particle-wall interaction reads 

Sd ° Ty 1 ojKo Tt ° Sd(v x , v y ), v x >—d (11) 
S d ° T T X o M^ 1 o T T ° S d (v x ,v y ), v x <-d, (12) 

where shifts of different sign are used for the upper (lower) wall to let the walls move into opposite directions. Other 
prescriptions to impose a shear will be investigated in more detail in the following section. 

In the simulations we set d = ±0.05, T = T u = T d = 1.0, N = 100 and p = 0.1. As we expected, we find a NSS 
with a linear shear profile along the x-direction (Fig. |J), where the drift velocity of the wall u w (*) is defined as the 
average between the in- and outgoing tangential velocities. The temperature profile is shown in Fig. ||, with the wall 
temperatures T w (+) defined as above. As can be seen in the plots, none of these values correspond to the parameters 
T (*) or d. Nevertheless, we obtain a linear shear profile u x {y) and an almost quadratic temperature profile T(y), as 
predicted by hydrodynamics [ 15[ . 

For a comparison of experimental and theoretical viscosity we follow the same procedure as above, i.e., we estimate 
the experimental shear rate du x (y)/dy by a linear least square fit u x (y) = jy, again discarding the outermost layers. 
Through the measured momentum transfer from wall to wall n the experimental viscosity is given as 

f] exp = n/7, (13) 

and the theoretical value as calculated by the Enskog theory has the form p3| , p^ | 



(Model I) 

( v ' X i v 'y) = 
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1 



m = 1.022: 



6n + 0.8729(6n) 2 x 



(14) 



Again b = n/2 denotes the second virial coefficient, \ 1S given by Eq.(^) and we use the arithmetic mean of the viscosity 
over layers 3 - 18 to compute the theoretical viscosity, i.e., r)th = ^-/^layersYl r ll- Table || shows good agreement 
between the experimental and the theoretical values for different numbers of particles, so again the thermostating 
mechanism leads to the correct macroscopic behavior. However, the discontinuities in the i^- velocity distribution of 
the scattered particles should be noticed (see FigJ^). This time these discontinuities do not diminish or disappear in 
the hydrodynamic limit. 



IV. ENTROPY PRODUCTION AND PHASE-SPACE CONTRACTION 



Having a deterministic and reversible dynamics at hand we can now turn to properties beyond the usual hy- 
drodynamic ones and, in particular, investigate the conjectured identity between phase space contraction rate and 



thermodynamic entropy production in the light of our formalism [ p3| , [24[]ig|jl5| , |l6| ,25|-|28|]. In an isolated macroscopic 
system the entropy is a thermodynamic potential and therefore plays the central role in determining the time evolution 
and the final equilibrium state. Yet, its microscopic interpretation out of equilibrium remains controversial (see, e.g., 
J|5],[l6]]) and the situation is even much less clear in NSS. For the class of models where thermostating is ensured by 
friction coefficients an exact equality between entropy production and phase space contraction rate in NSS has been 
inferred on the basis of a global balance between the system and the reservoir [p3|J^j9|,pl| . For the Chernov- Lebowitz 
model an approximate equality has also been found fl^ , |l6|| . Still, it is not clear under which circumstances this 
relation holds in general |l^,|l^,^ 53 . 



A. Equilibrium State 



We begin with the simplest case of equilibrium described in section [□], where the thermodynamic entropy production 
R eq vanishes. The bulk dynamics being Hamiltonian, phase-space contraction can only occur during collisions with 
a wall. Since these collisions take place ' instantaneously', we ignore the bulk particles and restrict ourselves to the 
compression related to a single collision during the time interval dt. The phase space contraction is then given by 
the ratio of the one-particle phase-space volume after the collision (dx' dy' dv' x dv' y ) to the one before the collision, 
(dxdydv x dv y ), and can thus be obtained from the Jacobi determinant of the scattering process. One easily sees that 
dx' = dx and \dy'/dy\ = \v' y dt/v y dt\ |l^JT^]. Furthermore, 



dv' x dv' y 



dv x dv y 



dT dM^V dT 



dv x dv y dxdy 



dv' x dv' y 



dT 



dv x dv v 



dT 



dv'dv' 



exp 



V' x 2 +V'y 2 



2T 



(15) 



where step two follows from the phase-space conservation of A4/A4 1 , and the last line is obtained from Eqs. 
Hence, in a particle-wall collision the phase-space volume is changed by a factor of 



dv' x dv' y dx' dy' 




= exp 1 


dv x dvydxdy 



v' 2 — v 2 — v 2 

y x u y 



2T 



(16) 



The mean exponential rate of compression of the phase space volume per unit time is thus given by 

dv' x dv' y dx' dy' 



P 



< In 



dv x dvydxdy 



>=((v 2 x + v 2 y -v' 2 -v' 2 )/2T) 



(17) 



where the brackets <> denote a time average over all collisions at the top and the bottom walls. In equilibrium the 
in- and outgoing fluxes associated to these collisions have the same statistical properties, so P eq sums up to zero and 



P, 



Req 1 



(18) 



This is fully confirmed by the simulations. 
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B. NSS 



The thermodynamic entropy production a per unit volume of our system in NSS is given by the Onsager form [[15|]S|] 

, , II du x T . . d ( 1 \ , „. 



T dy dy \T 

where II is the x-momentum flux in the negative y-direction, and J(y) is the heat flux in the positive y-direction 

1. Heat Flow 



Imposing only a temperature gradient on our system like in section [II A , the first term in Eq. (|T^) is identical to 
zero and the total entropy production R in the steady state is then 



R = o-dr= J/Tds = JUTl + J«/T« (20) 

J Volume J Surf ace 

The right hand side of Eq.(pO|) is the outward entropy flux J w /T w across the walls of the container. Note that there 
is no temperature slip at the walls with respect to the correctly defined temperature values, as indicated by the 
simulation results in Fig s.0 and |. 

On the other hand, the exponential phase-space contraction rate now reads 

P = +4~ v * ~ v 'y)l 2T ^) u + + v l - < 2 - v 'y)/ 2T d) d > (21) 

where we averaged over the upper and the lower wall separately. Since in NSS J™ — ((v 2 + v 2 . — v' 2 — v' y 2 )/2^ = 

~Jw = — (( v x + v y ~ v x ~ v y)/ 2 )d> ^ e rat i° s 01 entropy production to exponential phase space contraction rate 
reduce to 

Tf /d T u / d 

(22) 



— u/d rpu/d 



In the hydrodynamic limit the in- and the outgoing fluxes approach local equilibrium, implying T^ d ~ ~ r p^ l / d 



rpu/d £ or ^q^Jj wa u s Therefore, the ratios in Eq.(|22]) should go to unity. The numerical results in Table III confirm 
this expectation, leading to a good agreement between entropy production and exponential phase-space contraction 
rate. 



2. Shear Flow 

We follow the same procedure as in the preceding section. For a stationary shear flow the hydrodynamic entropy 



production a per unit volume in Eq. (19) can be written as |15|]8|| 

*>-fS + '<">s(f)- n s(T)- < 23 > 



dy dy \T J dy 

The second step in Eq.(|23|) follows from the fact that in NSS XdT/dy = J(y) = Hu x (y). The total entropy production 
R in the steady shear flow state is then 115|,n6l 



R= I adr= Wu/Tds 

J Volume J Surface 

J w /T w = 2L 2 II (u w /L) /T w = L 2 n 7 /T w . (24) 



In the macroscopic formulation of irreversible thermodynamics Eq. (|24|) is interpreted as an equality, in the stationary 
state, between the entropy produced in the interior and the entropy flow carried across the walls. Our shift map Sd 
mimics moving walls with drift velocities ±u w . The work performed at these walls is converted by the viscous bulk 
into heat and then again absorbed by the walls which now act as infinite thermal reservoirs. By imagining that the 
walls act as an 'equilibrium' thermal bath at temperature T w , R can be interpreted as their entropy increase rate. 



G 



For Model I (Eqs.(|ll])([l2|)) the mean exponential phase-space contraction rate takes the form (dSd/ dv x dv y = 1) 

P=-(K 2 + v' 2 -vl- v 2 y - 2d(v' x + v x )]/2T) , (25) 
whereas the entropy production is given by 

J w /T w = - (K 2 + v' 2 -v 2 x - v 2 - <v' x > 2 + <v x > 2 }/2T w ) . (26) 



In Table IV the ratios of L 2 H^/J W as obtained from the simulations are reported and the relation of phase-space 
contraction rate to entropy production is subsequently checked. Whereas the equality between entropy production 
and entropy flow via heat transfer is confirmed, we observe a significant difference between entropy production and 
phase-space contraction which subsists in the hydrodynamic limit. This mismatch is perhaps not so unexpected 
in view of the distorted outgoing fluxes (see Fig.^). Nevertheless, one could argue that the fine structure of these 
distributions may depend on the specific characteristics of the baker map chosen to model the collision process. We 
therefore also considered the standard map 



with the parameter k = 100 to ensure that we are in the hyperbolic regime f ff0|]41| ]. We found that the discrepancy 
between entropy production and phase-space contraction rate remains: although the deterministic map now seems to 
be chaotic enough to smooth out the fine structure of the outgoing densities, the discontinuity at d survives. Actually, 
as long as Model I is adopted it becomes clear that in NSS there will always be more out- than ingoing particles with 
v x > d at the upper wall (and with v x < d at the lower wall). Thus, the Gaussian halves in Figj^ (b) will never match 
to a full Gaussian even in the hydrodynamic limit, and $i and <j> will never come close to a local equilibrium. To 
circumvent this problem we modify the map T and investigate the following model: 



(Model II) 



erf 



\v x \Td)/V2T ±erf(d/V2T) 

T ± (v x ,v v ) = I L - j -= ,exp(-v 2 /2T) I (28) 

1 l±erf(d/V2T) " ' 

(v' x ,v' y ) =T~ 1 oMoT-{v x ,v y ), v x >0 (29) 
{v' x ,v' y ) = TZ 1 o M- 1 o T +(v x ,vy), v x < 0. 

This model is also time-reversible, but in contrast to the former one no particle changes its tangential direction during 
the scattering. There is still a gap in the outgoing distribution of Fig.(|g) (b), however, simulations show that this gap 
disappears in the hydrodynamic limit thus bringing the in- and the outgoing distributions close to local equilibrium. 
Furthermore, we note that whereas we were not able to give a relation between the parameter d and the actual wall 
velocity u w for Model I, in case of Model II u w converges to d in the hydrodynamic limit. For this reason we chose 
d = 0.5 in the following, since this value yields the same order of the wall velocity as d = 0.1 for Model I. Proceeding 
now to the phase-space contraction rate we find that it takes the form 



^ _ m„ 1 ± erf(d/v/2T) 
1 =f erf(d/v / 2T) 

where n± are the collision rates for positive and negative tangential velocities, and the additional term (cf. Eq.([25|)) 
results from the different denominators in Eq.(2p). Note that one has to average over the upper and the lower wall 
separately. Again we compare R and P (Table V[), but although the outgoing flux approaches now a Gaussian in the 
hydrodynamic limit the two quantities still do not match. This result can be understood in more detail by rearranging 
the terms in Eq. j30|) as 

P U ' d = ~ (K 2 + v' 2 -v 2 x - v 2 y - <v' x > 2 + <v x > 2 ]/2T) u/d 

- ([< v' x > 2 - < v x > 2 -2d(v' x + v x )]/2T) u/d - (n + - n_) In . (31) 



Since T — > T w in the hydrodynamic limit, the first term clearly corresponds to the entropy production Eq.(|26j). 
However, the second and the third terms provide additional contributions. For u w — > d and d — » they are both 
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of order d 2 and can be interpreted as a phase space contraction due to a friction parallel to the walls |3(J. These 
two terms apparently depend on the specific modeling of the collision process at the wall. They may physically be 
interpreted as representing certain properties of a wall, like a roughness, or an anisotropy. Actually, the second term 
already appeared in Model I, see Eq. (|25|) . The price we had to pay in Model II for the fluxes getting close to local 
equilibrium is the additional third term in Eq.(|3l|), which does not compensate the second one. 

The foregoing analysis shows clearly what to do to get rid of the additional term in Eq.(^o|): we have to use the 
same forward and backward transformations T ± in Eqs.( p8|]29| ). If one still wants to transform onto a full Gaussian 
in the hydrodynamic limit time-reversibility has to be given up. This leads us to propose the model 

(Model III) 

(v' x ,v' y ) = T- 1 oMoT*(v x ,v y ). (33) 

which is still deterministic, but no longer time-reversible. The phase-space contraction is now given as 

P = {Wx + v'y ~v 2 x -v 2 y - 2d(v' x - v x )]/2T) . (34) 

Fig. ^ shows that the in- and the outgoing fluxes are getting close to local equilibrium, implying that the velocity 
of the wall u w goes to d and the wall temperature T w goes to T. Consequently, Eq.(|34|) should converge to the 
correct thermodynamic entropy production of Eq.(^6|) in the hydrodynamic limit, and this is indeed what we observe 
in Table This implies that time-reversibility does not appear to be an essential ingredient for having a relation 
between phase space contraction and entropy production, as was already stated in Refs. (ljjl^,^]. We remark that 
we consider the lack of time-reversibility in Model III rather as a technical difficulty of how we define our scattering 
rules than a fundamental property of this model. 



erf (v x - d)/V2T 



,exp(-^/2T) 



(32) 



V. CONCLUSION 



We have applied a novel thermostating mechanism to an interacting many-particle system. Under this formalism 
the system is thermalized through scattering at the boundaries while the bulk is left Hamiltonian. We have shown 
how this deterministic and time-reversible thermostating mechanism is related to conventional stochastic boundary 
conditions. For a two-dimensional system of hard disks, this thermostat yields a stationary nonequilibrium heat 
or shear flow state. Transport coefficients obtained from computer simulations, such as thermal conductivity and 
viscosity, agree with the values obtained from Enskog's theory. 

Having a time-reversible and deterministic system we also examined the relation between microscopic reversibility 
and macroscopic irreversibility in terms of entropy production. We find that entropy production and exponential 
phase-space contraction rate do in general not agree. When the NSS is created by a temperature gradient both 
quantities converge in the hydrodynamic limit. By subjecting the system to a shear we examined three different 
versions of scattering rules, of which one (Model III) produced an agreement. 

Our results indicate that neither time-reversibility nor the existence of a local thermodynamic equilibrium at the 
walls are sufficient conditions for obtaining an identity between phase space contraction and entropy production. A 
class of systems where such an identity is guaranteed by default are the ones thermostated by velocity-dependent 
friction coefficients |^J^,|ll|. We suggest that in general, that is, by using other ways of deterministic and time- 
reversible thermostating, such an identity may not necessarily exist. We would expect the same to hold for any 
system where the interaction between bulk and reservoir depends on the details of the microscopic scattering rules. 

As a next step it would be important to compute the spectrum of Lyapunov exponents for the models presented 
in this paper. This would enable to check, for example, the validity of formulas which express transport coefficients 
in terms of sums of Lyapunov exponents |l9|j2C| ], and the existence of a so-called conjugate pairing rule of Lyapunov 
exponents Moreover, it would be interesting to verify the fluctuation theorem, as it has been done recently 

for the Chernov-Lebowitz model ]37| ]. 
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FIG. 1. Temperature (-) and variance (- -)/(-•) profiles, T u — 2, T d = 1 (*). (+) denotes the wall temperatures T w as defined 
in the text. 




0.115 0.12 0.125 0.13 0.135 0.14 0.145 0.15 

number density 

FIG. 2. The profile of the number density, T u = 2, T d = 1 . 
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FIG. 5. Velocity variance in the a;-direction (- .) and in the y-direction (- -), temperature (-), measured wall temperature 
T w (+) and 'parametrical' temperature T (*). 
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FIG. 7. Velocity distributions of the in- and outgoing particles at the upper wall for the shear flow case (Model I, standard 
map): a) v?, b) vT\ c) <\ d) v° ut . 
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TABLE I. Comparison of theoretical and experimental heat conductivity \ £xp /\th 

N=100 N=200 N=400 N=800 

AT=1.5-1 0.904 1.009 1.003 1.062 

AT=2.0-1 0.887 0.950 1.021 1.051 



TABLE II. Comparison of theoretical and experimental viscosity rj exp /rith 

N=100 N=200 N=400 N=800 

d=0.05 0.9616 0.9904 1.0081 1.0382 

d=0.1 0.9702 1.001 1.0226 1.0232 



TABLE III. Comparison of entropy production and exp. phase-space contraction rate R/P, heat flow 







N=100 


N=200 


N=400 


N=800 


T u = 


1.5 


1.0814 


1.0762 


1.0614 


1.0508 


T d = 


1 


0.8948 


0.9170 


0.9273 


0.9439 


T u = 


2 


1.1313 


1.1110 


1.0985 


1.0765 


T d = 


1 


0.8122 


0.8412 


0.8633 


0.8886 



TABLE IV. Comparison of entropy production and exp. phase-space contraction rate R/P, shear flow, Model 

I 







N=100 


N=200 


N=400 


N=800 


rf=0.05 
baker map 


L A W 1 /J W 
R/P 


0.9816 
0.6761 


0.9669 
0.5882 


0.9765 
0.5023 


0.9962 
0.4230 


d=0.1 
baker map 


L z U~t/J w 
R/P 


0.9829 
0.6457 


0.9664 
0.5761 


0.9497 
0.4934 


0.9588 
0.4275 


d=0.1 
standard map 


L'Uy/J w 
R/P 


1.0255 
0.4622 


1.0424 
0.3873 


1.0435 
0.3008 


1.0833 
0.2417 
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TABLE V. Comparison of entropy production and exp. phase-space contraction rate R/P, shear flow, Model 



II and III 







N=100 


N=200 


N=400 


N=800 


Model II 




0.9842 


1.0154 


1.0002 


1.0255 


d=0.5 


R/P 


0.1858 


0.1682 


0.1398 


0.1127 


Model III 


z/n 7 /j ro 


1.0164 


1.0046 


1.0005 


1.0039 


d=0.5 


R/P 


0.8452 


0.8785 


0.9051 


0.9619 
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